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Surface plasmons in graphene offer a compelling route to many useful pho¬ 
tonic technologies [1—3]. As a plasmonic material, graphene offers several in¬ 
triguing properties, such as excellent electro-optic tunability[4], crystalline sta¬ 
bility, large optical nonlinearities [5] and extremely high electromagnetic field 
concentration[6]. As such, recent demonstrations of surface plasmon excitation 
in graphene using near-field scattering of infrared light [7, 8] have received intense 
interest. Here we present an all-optical plasmon coupling scheme which takes 
advantage of the intrinsic nonlinear optical response of graphene. Free-space, 
visible light pulses are used to generate surface plasmons in a planar graphene 
sheet using difference frequency wave mixing to match both the wavevector and 
energy of the surface wave. By carefully controlling the phase-matching condi¬ 
tions, we show that one can excite surface plasmons with a defined wavevector 
and direction across a large frequency range, with an estimated photon efficiency 
in our experiments approaching 1CT 5 . 

Graphene has attracted significant interest in recent years as a unique optical material. In 
particular, it has been predicted and experimentally shown that graphene can support very 
highly confined surface plasmon modes[l, 9], with an electrically tunable dispersion[7, 8]. 
Despite these promising discoveries, the burgeoning field of graphene plasmonics has some 
serious obstacles to overcome if it is to progress from the proof-of-principle stage. Problems 
arise due to the small wavelength of the surface plasmons, two orders of magnitude smaller 
than light of the same frequency. This has led to the development of specialised measurement 
techniques, most of which use infrared light and geometries with scattering resonances [10— 
12] or near-field sources [7, 8] to excite graphene surface plasmons. However, the far-infrared 
regime, in which graphene plasmons are predicted to have long lifetimes, lacks developed 
sources and detectors compared to the visible regime. Alternative approaches, such as the 
manipulation of surface acoustic waves to couple to the graphene surface plasmons [13, 14], 
therefore hold promise. 

Particularity desirable is the potential to excite a plasmon eigenstate with a singular 
energy, momentum and direction, vital for many future applications, including plasmonic 
circuits. In this respect, very recent progress has been made, with the development of 
carefully designed nano antennas which can locally excite and direct surface plasmons in 
graphene [11], Here, the combination of infrared source frequency and nanoantenna dimen- 
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sions determine the frequency, wavevector and direction of the surface plasmons generated. 
In this article, we investigate a competing approach that embodies many of these desirable 
aspects of directivity without requiring careful nanofabrication of antennas. This all-optical 
approach can access a distinctly broad frequency range, even down to the far infrared. We 
coherently excite surface plasmons using two visible frequency free-space beams via differ¬ 
ence frequency generation (DFG), an effect which we monitor through changes in reflectance, 
and can tune the frequency and wavevector of the surface plasmon through careful adjust¬ 
ment of incident light sources. This potential to excite and detect plasmons purely with 
free space optics, and at frequencies different than that of the plasmons themselves, has the 
potential to significantly expand the technological possibilities for graphene plasmonics. 

The intrinsic nonlinear interactions of graphene with light are surprisingly large[5, 15- 



FIG. 1: The nonlinear coupling scheme illustrated on a dispersion diagram. The DFG of the pump 
(green arrow) and probe (orange arrow) allows access to wavevectors outside of the light line (red 
line). This permits phase-matching to the surface plasmon modes in graphene (blue line). The pink 
line illustrates a region that can be interrogated by altering the pump wavelength from 615 nm to 
545 nm with the probe wavelength fixed at 615 nm. (Inset) The experimental arrangement used 
to excite surface plasmons on graphene. 
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18]. Moreover, large enhancements of nonlinear optical effects are predicted by the presence 
of highly confined plasmons in graphene[19, 20]. It seems intuitive, then, to attempt the 
converse: to use the nonlinear interaction between optical fields to resonantly drive surface 
plasmons. This kind of approach has been demonstrated for thin metallic films [21, 22], and 
has been recently proposed for graphene, with various coupling schemes for the difference 
frequency mixing of infra-red light in graphene clad waveguide structures suggested [23, 24], 
Similar in concept, figure 1 shows our nonlinear coupling scheme illustrated on a dispersion 
diagram. By illuminating the graphene with two intense laser pulses with well-defined angles 
of incidence but different frequency, labeled here f pum p and f pro be, one can phase match 
both the frequency and wavevector, k , of the surface plasmon. This wave mixing process 
is a second order nonlinear effect, normally forbidden in centro-symmetric crystals[25], but 
possible in graphene because of the distinctively non-local, spatial character of the interaction 
[20]. The inset in fig. 1 shows the experimental arrangement used. An identical pair of 
optical parametric amplifiers (OPAs), pumped by an amplified femtosecond laser system, 
generate the 100 fs-pulses at a repetition rate of 1 kHz. The wavelengths of the two OPAs 
are selected independently, and the beams are directed to the sample. The experiments 
are carried out in a non-collinear geometry, using two beams incident on the samples at 
angles 9 pump and 9 pro b e • Oblique angles provide sufficient in-plane momentum to match to 
the surface plasmon, as illustrated in fig. 1. The incident beams are weakly focused on 
the sample using 30 cm focal length lenses, giving rise to a very small uncertainty in angle 
~ 0.017 rad, and a similarly negligible uncertainty for the in-plane wavevectors. Sets of 
half-waveplates and polarizers determine both the average power and polarization, with 
the polarization set such that the electric vector of the light is in the plane of incidence 
(transverse magnetic polarized). The pump pulse fluence, <f>, used is typically in the range 
<f> ~ 0.1 — 0.2 mJ/cm 2 , with a pump spot size on the sample of ~ 300 //m radius. This pump 
fluence is an order of magnitude less than the photo-modification threshold for graphene [26], 
and the probe fluence is typically two orders of magnitude smaller still. We measure the 
differential reflection of the probe beam defined as A R/R — (R — Rq)/Ro where R and Ro 
are the reflections with and without the presence of the pump pulse, respectively. In order 
to isolate the nonlinear reflection signal, we vary the temporal overlap of the two pulses 
using a motorized delay stage. 

For optical excitation pulses, one expects optical nonlinearity arising due to saturable 
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absorption caused by Pauli blocking of interband transitions[27]. A typical measurement 
of the temporal dynamics recorded for this process (A pump = 547 nm, \ pro b e = 615 nm) 
is shown by the black curve in fig. 2. The asymmetric line shape of the signal is due to 
the relaxation dynamics of the excited electrons cooling[28, 29], with significant temporal 
broadening caused by the spatial overlap of the non-collinear beam spots. Note that there 
is no appreciable signal from the quartz substrate (See Supplementary Information (SI) 

fig. SI). 

For non-degenerate pump and probe beams, in addition to (incoherent) saturable absorp¬ 
tion effects, one can expect (coherent) wave mixing signals. This coherent contribution to 
the probe reflection is expected to be significantly enhanced when the difference frequency 
held generated by the pump and probe matches that of the graphene surface plasmons. 
This is analogous to that of a stimulated Raman process, corresponding to a transfer of 
energy from pump to probe pulses[25] via the generation of surface plasmons. An example 
of the recorded temporal dynamics under such a resonant condition is presented in fig. 2. 
Comparing the two curves in this figure, we see that the “non resonant” signal (i.e. when 
one is not phase matching to plasmon excitation) gives rise to an asymmetric lineshape 
representative of carrier cooling dynamics. Under “resonant” conditions (i.e. when phase 
matching conditions are satisfied) we observe a fast additional contribution to the signal, 
giving rise to a more symmetric lineshape, as one would expect for a coherent signal. For 
certain experimental geometries and excitation huences, signal enhancements of up to ap¬ 
proximately x4 can be observed (see supplementary figure S3). It should be noted that, 
depending on efficiencies, it may be possible to isolate the coherent signal using a heterodyne 
detection scheme[30], which could also allow detection of a plasmon in a different spatial 
position than generated. 

In order to observe the presence of the coherent signal, we vary the difference frequency 
in order to isolate any resonant, coherent conditions. In experiment, the pump wavelength is 
varied from 615 nm to 545 nm, with the probe wavelength set at 615 nm, allowing difference 
frequencies ranging from 0 to 60 THz. In this way, it is possible to interrogate a section 
of the surface plasmon dispersion, for example the region illustrated by the pink line in 
fig. 1. Note that we normalize the signal by pump fluence in order to remove artefacts due 
to power variation[28]. By altering the experimental geometry, we investigate here three 
different regions of the dispersion diagram corresponding to (0 pump = 55 °,0 probe = 45°), 
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Figure 3 shows the results of these three measurement geometries, superimposed on the 
surface plasmon dispersion (black line). The dispersion was calculated according to the 
model outlined in ref. [31], with the SiC >2 substrate phonon frequencies as given and a Fermi 
energy of Ef = 0.5 eV. This Fermi energy is larger than the expected intrinsic doping of 
our graphene samples (see Methods for sample details), which we attribute to a significantly 
raised electron temperature expected under illumination by ultrafast pulses (see SI, fig. S2). 
Hybridization with the substrate phonons leads to four branches[31, 32], The overlaid colour 



delay time (ps) 


FIG. 2: Differential reflection normalized to average pump fluence as a function of temporal overlap 
for the geometry 0 pU mp = 15°, 0 pro b e = 125°. At zero delay time, both the pump and probe pulses 
arrive simultaneously, leading to a nonlinear change in the probe reflection. Two curves are shown: 
The black curve labelled “non-resonant” shows a typical time asymmetric mesurement when the 
difference frequency produced by the pump and probe (61.2 THz) does not coincide with a surface 
plasmon energy state. The red curve shows an additional fast symmetric contribution to the 
recorded reflection signal when the differency frequency matches the energy of a graphene surface 
plasmon (23.8 THz). 
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plots are placed on the diagram so that the maximum differential reflection signal achieved in 
each delay-scan corresponds to the difference frequency and wavevector of the dataset. The 
grey shading around the plasmon dispersion curve indicates the expected spectral broadening 
of the signals (~ 7.5 THz) due to the finite bandwidth of ~ 100 fs-pulses. 

Near the regions defined by the surface plasmon dispersion in graphene, we observe 
clear enhancement in the differential reflection. The assignment of the spectral features to 
surface plasmon excitation is further supported by the polarization dependence of the signal 



FIG. 3: Plots of differential reflection normalized to pump fluence for the three different experimen¬ 
tal geometries, superimposed on the graphene surface plasmon-phonon dispersion, (black lines). 
Angles used are: (a) 9p U mp = 55 , dprobe = 45 , (b) Opump = 50 , 9probe = 70 and (c) 9p U mp = 15 , 
9probe = 125°. The intraband transition threshold and light line (dotted lines) are labelled on the 
diagram. The grey regions indicate the spectral resolution of our experiment. 
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(see fig. S5). The observation of these resonant features over the incoherent background is 
strongly dependent on the magnitudes of both pump and probe intensities (see SI figures 
S2 and S3). For larger difference frequencies, up to 150 THz, we do not observe any further 
resonance features in our spectra (see fig. S4). 

The lower branch of the plasmon dispersion relation gives rise to the largest mixing 
signals for the low wavevector phase-matching (angles (a) and (b) in fig. 3), while the 
upper frequency branches give rise to largest signals for the high wavevector (angle (c) in 
fig. 3) region. Whilst we observe clear resonance features in all three of these experimental 
geometries, we also observe a change in sign of the signal between low (figure 3(a) and (b)) 
and high (figure 3(c)) wavevector regions. The absolute differential reflectivity signal size 
also increases with increasing wavevector. 

To understand the origin of these coupling behaviours, we have developed a simple theo¬ 
retical model that captures the salient features of this nonlinear reflection and generation of 
plasmons, and we briefly summarize the ideas here (details of the model are presented in the 
supplementary material). In general, equations for the electromagnetic boundary conditions 
at the air-graphene-substrate interface relate the wavevector and frequency dependent re¬ 
flection and transmission coefficients r(k,u>),t(k,u )) to the graphene current density J(k,u>). 
The current density, on the other hand, can be written in terms of the electric field via 
conductivity functions, which allows the equations to be solved in terms of fields alone. 
Nonlinear contributions imply that J(k,oj) depends on fields at other wavevectors and fre¬ 
quencies, which couple the various reflection and transmission coefficients together. For a 
second-order conductivity cb 2 ), we find that the probe transmission depends on the pump 
via 

f {L) 

+ _ probe 

I 2 / 

pump | -bpump 

with an analogous equation for the pump transmission (expressions for r are more involved 
but are directly related to t, see supplementary materials). Here t^ obe is the linear trans¬ 
mission coefficient, A is a function of linear optical properties and beam angles, and for 
notational simplicity, dependencies on k,u are implicit. In Fig. 4, we plot the numerical 
solution for the differential probe reflectance, normalized by fluence, for the simplified case 
of continuous plane-wave pump and probe beams. While this simple model ignores the 
non-equilibrium nature of the excitation, as observed in experiment, we show below that 
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it is sufficient to describe some of the salient features of our results. Similar to Fig. 3, the 
differential reflectance is plotted versus difference frequency and in-plane wavevector, whose 
values are scanned by continuously varying the probe incidence angle and pump wavelength. 
The probe wavelength and pump angles are fixed at 615 nm and 50°, respectively, and the 
pump and probe intensities are chosen to be 10 and 0.1 W//mi 2 , to closely correspond to 
the configuration in Fig. 3(b). It can be seen that the simulation qualitatively produces the 
main features of Fig. 3. In particular, the change in the sign of differential reflectance at the 
Brewster angle is clearly observed, as is the enhancement of the signal when the difference 
frequency and wavevector align with the plasmon dispersion relation. 

The model also reproduces some of the main features arising from different coupling 
efficiencies to different bands (fig. 4). Generally, the highest coupling efficiency occurs for 
the dispersion regions which are most ‘plasmon-likc’ in origin (see additionally fig. S8). 
This is most obvious comparing the data in fig. 3(c) and fig. 4(c), where the coupling to the 
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FIG. 4: The numerical solution for the differential probe reflectance normalized by fluence, calcu¬ 
lated using the model outlined in the SI. The white dotted lines indicate the region of the dispersion 
relation probed by the experimental geometries shown in fig. 3. 
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upper bands is much stronger than the lowest band in both model and experiment. For lower 
wavevector cases, the coupling to the highest band is overestimated in the model compared 
to the experiment. This could possibly be caused by frequency-dependant losses in the 
graphene sheet unaccounted for in our simple model. The model additionally reproduces 
the increasing absolute signal strength with increasing wavevector observed in experiment, 
which is a consequence of both larger changes in the reflection coefficient for a corresponding 
change in absorption for higher angles, and due to spatial dispersion[20] in the signal. Indeed, 
it can be shown that the magnitude of enhancement is proportional to the square of the 
plasmon quality factor, Q 2 (see SI), in agreement with predictions from ref. [23]. 

In addition to the surface plasmon resonance conditions, for the highest wavevector region 
in fig. 3(c) there is an additional resonant enhancement found at low frequencies < 3 THz 
measured in experiment that is not reproduced in our model (fig. 4(c)). The position of 
this peak lies within the expected region of intraband transitions in graphene, indicated by 
the dotted line in fig. 3. This feature is also largely polarization independent, unlike the 
enhancements we attribute to surface plasmon coupling (see SI, fig. S5). 

In principle, Eq. 1 can be inverted to allow an experimental determination of the nonlin¬ 
ear conductivity a^ 2 \ given transmission or reflection data. This is difficult in the present 
setup, in part given the broad bandwidth of the pulses, uncertainty over some system 
parameters, and difficulty of investigating a large number of angles to quantify possible 
wavevector and frequency dependence of cR 2 b However, as an estimate, we take the sim¬ 
plest possible model, in which the effective nonlinear susceptibility is frequency and 
wavevector-independent. This corresponds to a nonlinear conductivity function obeying 
cr l ' 2 \u) = i^^ijjJprobe)^/ui pr obe)i where the value at the (fixed) probe frequency represents 
a single fitting parameter. We find that a value of \cr l ' 2 ' ) (ujp ro be)\ ~ 2.4 x 1CT 12 A ■ m/V 2 
produces the same peak signal as observed in Fig. 3(b). It should be emphasized that 
this represents a rather conservative estimate of a^ 2 \uj). In particular, the mobility of 
H « 2000 cm 2 /Vs corresponds to a plasmon linewidth of 7 = 27 T x 1.6 THz that is nar¬ 
rower than the measurement bandwidth, indicating that only a fraction of the pulse can 
efficiently excite plasmons. Reducing measurement bandwidth could therefore give rise to 
greater coupling to the surface plasmons, while also reducing the effects of non-equilibrium 
carriers on the measurements. While a comparison to a bulk nonlinear crystal is not di¬ 
rectly meaningful, it is nonetheless interesting to note that a bulk nonlinear crystal with 
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the thickness t ~ 0.3 nm of a graphene layer would require a nonlinear susceptibility of 
^( 2 ) rsj (hJprobe) I /(^o^probe^) ~ 3 x 10 -7 m/V to produce the equivalent in-plane nonlinear 

currents. This value is approximately 3 orders of magnitude larger than in GaAs. 

Finally, from the inferred value of a^ 2 ' 1 and the input beam parameters, our model enables 
us to estimate the conversion efficiency rj of pump photons to plasmons (see SI). We find a 
value of r] « 6 x 1CT 6 , while noting that the actual conversion could be significantly higher 
with narrow pulses, again as the estimated value of does not account for the large pulse 
bandwidth. We note that this experimentally obtained value of r/ is of the same order as 
predicted in ref. [23], once adjusted for our experimental parameters. 

In summary, surface plasmon generation in a planar single layer graphene crystal has 
been demonstrated using nonlinear wave mixing of visible frequency light. We observe 
enhancements in the nonlinear, differential reflection of light for various frequencies and 
phase-matching conditions which agree well with the surface plasmon-phonon dispersion for 
graphene on quartz substrates. In contrast to near-held infrared scattering[7, 8], our ap¬ 
proach may be used to access an distinctly broad frequency range, including the ordinarily 
hard to reach THz “gap” [33]. Moreover, by careful manipulating the phase-matching condi¬ 
tions, we show that one can generate surface plasmons with a defined wavevector, with an 
efficiency approaching 1CT 5 . This efficiency by no means represents a fundamental limit, and 
we believe that it could in principle be pushed towards a 10~ 2 level with future adjustments, 
such as increasing the surface plasmon Q factor from ~ 5 to ~ 30 with lattice-matched 
hBN substrates[34], equalizing the intensities of the pump and probe beams (see SI, fig. 
S3), or the use of narrower bandwidth pulses. Moreover, in principle, our approach could 
be extended to higher or lower frequencies, regions that are generally hard to access using 
current approaches [2]. 

METHODS 

Sample Preparation 

Samples for our experiments are fabricated from commercially grown CVD graphene on 
copper foil (graphene supermarket). Transfer to quartz substrates was performed in house 
via a standard metal etching and float technique using Ammonium persulfate to etch the cop- 
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per and PMMA as a support structure. Combined resistance and Raman spectroscopy[35] 
give an estimated mobility of the samples of around 2000 cm 2 /Vs and a natural Fermi energy 
of ~ 300 meV. Raman imaging indicates that the graphene is nominally single layer, with 
~ 80% coverage of the substrate. 
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SUBSTRATE RESPONSE 

Differential reflection was recorded as a function of delay time for both the graphene 
on quartz and for the bare quartz substrate to ascertain any contribution to the nonlinear 
response from the substrate. The experimental parameters for these data were set to the 
resonant condition of figure 3(b), with 9 pump = 50°, 0 pro b e = 70° and the difference frequency 
set to 10 THz. The measurements are compared in figure SI, and negligible signal, compared 
to the resonant measurement from graphene, is observed. 



FIG. SI: Differential reflection normalized to fluence as a function of temporal overlap for: 
Graphene on quartz (black) and the bare quartz substrate (red) at d pum p = 50°, 6 pro b e = 70°. 
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INCREASING PUMP FLUENCE 


For pulsed fluence ~ 0.1 mJ/cm 2 we expect to generate an electron temperature of 
~ 1000 K[SI]. This means that we are probing a very non-equilibrium electron distribution. 
In order to investigate the effect of this electron heating, we increased the pump fluence 
used in the experiment to ~ 1.1 mJ/cm 2 . These results are shown in fig. S2. We find 
that a higher fluence significantly suppresses the surface plasmon resonance features with 
respect to the background, off-resonance signal. We believe two factors contribute to this 
effect: firstly, an increased electron temperature will increase saturable absorption, the effect 
primarily responsible for the incoherent background signal. Secondly, due to the negative 
photoconductivity usually exhibited by graphene for pulsed femtosecond excitation[S2], one 
can expect increased losses and quenching of the surface plasmon, leading to broadening of 
the spectral features associated with their excitation. 

Evidence that the increased electron temperature also raises the effective Fermi energy of 
the sample can be inferred by comparing fig. S2(a) and fig. S2(b), where the resonant regions 
of differential reflectivity shift to higher frequencies in the high-fluence case, as would be 
expected for the graphene surface plasmon dispersion for a higher doping level. 
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FIG. S2: Differential reflection normalized to fluence as a function of temporal overlap for: low 
(left) and high (right) pump fluence at 9 purnp = 15° and 9 pro be = 125°. 
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REDUCING PUMP FLUENCE 


Since electron heating clearly has a negative impact in the experiment, one wants to min¬ 
imize the pump intensity. However, reducing the light intensity also reduces our difference 
frequency coupling efficiency, since surface plasmon generation here is a nonlinear process. 
However, a possible route to better isolating coherent signals is to reduce the pump beam 
intensity and increase the probe beam intensity, illuminating the sample with similar Ali¬ 
enees for both beams. This reduces the Pauli blocking of the probe induced by the pump 
beam, decreasing the signal due to saturable absorption decreases, while maintaining the 
efficiency of the difference frequency mixing process. 

Figures S3 (a) and S3(b) show differential reflection normalized to pump fluence for two 
difference frequencies, 0 THz (A p Urnp = 615 nm) and 12 THz (\ pump = 600 nm), measured for 
the angles 9 pump = 50° and 9 pro i, e = 70°. In this geometry we expect a resonant enhancement 
for a difference frequency ~ 12 THz due to plasmon excitation and no enhancement for 
0 THz, as previously measured in fig. 3(b). In figure S3 (a) we show a typical measurement 


(a) (b) 



FIG. S3: Differential reflection normalized to pump fluence as a function of temporal overlap of the 
pulses for: (a) high-fluence pump and low-fluence probe, and (b) comparable fluence in both pump 
and probe. The black lines indicate the case for the difference frequency of 0 THz (no resonant 
plasmon coupling) and the red lines are for a difference frequency of 12 THz (resonant plasmon 
Coupling). 0p Urn p — o5 , Oprobe — d) , A probe — 615 nm. 
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for a high-power pump beam (0.26 mJ/cnr) and a low power-probe beam (0.0028 mJ/cm 2 ). 
In this case, when the difference frequency matches that of the surface plasmon, a resonant 
change to the reflectivity by a factor of 1.6 occurs. For equal pump and probe fluences (~ 
0.07 mJ/cm 2 ), as shown in fig. S3(b), we observe a significant suppression of the background, 
non-resonant signal. The enhancement then measured from a non-resonant condition to the 
resonant surface plasmon excitation is increased to a factor of 3.9. 

We also note that in the low-power pump case (fig. S3(b)), the lineshape of the recorded 
temporal dynamics is far more symmetric than in the high-power case (fig. S3 (a)), which 
clearly exhibits the typical asymmetric lineshape indicative of incoherent carrier cooling 
dynamics. 

HIGHER DIFFERENCE FREQUENCIES 



FIG. S4: Plot of differential reflection normalized to pump fluence as a function of temporal overlap 
for 9 pump = 15° and 9 pro b e = 125° at higher frequencies. The colorscale has been scaled to fig. 3(c) 
for ease of comparison. 

We have preformed several measurements for a wider range of difference frequencies. In 
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FIG. S5: Differential reflection normalized to pump fluence as a function of temporal overlap for 
transverse magnetic (a) and transverse electric (b) polarizations. The resonances for the predicted 
surface plasmon frequencies are clearly suppressed (labelled), while the intraband resonance is 
largely unaffected. 

figure S4, we present measurements taken for 9 pump = 15° and 9 pro b e = 125°, where the pump 
wavelength was varied from 540 nm to 475 nm, with the probe wavelength fixed at 615 nm. 
This gives a difference frequency range from 70 THz to 140 THz. For these larger difference 
frequencies, we observe no resonance features above 70 THz, indicating there is no coherent 
coupling to higher frequency modes. 


POLARIZATION DEPENDENCE 

The experiment shown in fig. 3(c) was repeated with the pump and probe both polarized 
with the electric vector parallel to the graphene surface (transverse electric, TE polarized). 
Under these circumstances we expect surface plasmon excitation to be suppressed. The 
results are shown in figure S5. When illuminating with TE polarized light, there is a clear 
decrease in the nonlinear enhancement to the reflectivity for the 23-27 THz peak and the 
~ 45 THz peak. The peak at ~ 0 THz is of the same order in both polarization cases, 
suggesting that the nonlinearity arising due to the intraband transitions are less sensitive 
to the polarization of the light. This is attributed to how carriers respond to the EM held 
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in the plane of the surface. Since, for the intraband response the conductivity functions are 
in-plane, they make no distinction as to the out of plane component (which distinguishes 
TE and TM polarization). 

While the higher-frequency resonances are strongly suppressed, they are not absent com¬ 
pletely: we attribute this to imperfect polarization of the beams and inhomogeneity in the 
graphene sample. 


THEORETICAL MODEL 


We describe a theoretical model of two continuous-wave, free-space beams of frequencies 
uq i2 (without loss of generality, assume that uq > lo 2 ) interacting with graphene via a 
difference frequency generation process. The convention in our calculations to define the 
held polarizations and beam angles is illustrated in Fig. S7. The beams are taken to be 
incident from air (refractive index n ~ 1). Important to modeling this experiment is the 
inclusion of a frequency-dependent and complex refractive index of the substrate at low 
frequencies, in order to capture the lattice vibrations in silica and the resulting surface 
optical phonons. To do this, we take a simple dielectric response model based upon three 
transverse optical (TO) phonon modes [S3], 

3 


n 2 (uj) = Cpo + y] 


r 2 

Jj^TOj 


y ^ TO,j - ^ ~ iuiTOj ' 


(SI) 


From Ref. [S3], the high-frequency dielectric constant is taken to be = 2.4, while the 
TO phonon frequencies and oscillator weights are u>to = 2n x (13.44, 23.75, 33.84) THz 
and / = (0.7514,0.1503,0.6011), respectively. The damping rates are taken to be q^o = 
2n x (0.80,1.27,1.27) THz. The resulting real and imaginary parts of the refractive index, 
plotted in Fig. S6, approximately correspond to experimentally measured values [S4], In 
practice, this refractive index function is only relevant for the substrate response at the low 
difference frequency of = oq — oq, while for the high frequencies uq^ the response is nearly 
frequency-independent, n ~ ^/e^. 

In general, one can obtain equations relating the reflection and transmission coefficients 
to each other by enforcing electromagnetic boundary conditions (continuity of the normal 
electric displacement and tangential electric held) at the graphene interface. The solution 
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for the transmission coefficient of field i (i = 1, 2) is readily found to be 

t _ 2 sin 6j _ (pjs/e 0 ) sec fa sin 0, 2 

1 rii sin 9i + sin (pi Eu (nf sin 9i + rii sin (pi) 1 

where n* = n(u>i) denotes the substrate refractive index at the held frequency, while the 
reflection coefficient is related by r* = 1 — t t sin esc 0*. Here Eu are the incident held 
amplitudes, and 9 1: and & are the angles of the helds on the vacuum and substrate sides, 
respectively. p is = p(oJi, k ix ) is the graphene surface charge density at frequency uii and 
in-plane wavevector ki X = (cu*/c) cos 9i . Note that the hrst term on the right-hand side of 
Eq. (S2) reproduces the standard Fresnel coefficient in the absence of a graphene layer (p is = 
0). The angle of the transmitted held is related to the incident by Snell’s Law, cos 9i = 
n l cos (pi. 

The surface charge density can be related to the current density J in the graphene layer 
via the continuity equation, which in the Fourier domain reads 

Ps(w, k x ) {kx/JxipJi kx)• (S3) 



FIG. S6: Real (black) and imaginary (red) parts of the refractive index of silica versus free-space 
wavelength (A = 2-kc/uj ), based upon the model of Eq. (SI). 
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At the same time, the current density can be related to the electric holds via conductivity 
functions. We are particularly interested in the case of difference frequency generation, 
where the held produced at the difference frequency and wavevector uj 3 — uq — u> 2 and 
k 3x = k\ x — k 2x is resonantly enhanced by aligning them with the plasmon dispersion of 
graphene co p (k x ). This motivates a truncated model in which we include only the linear 
and second-order conductivities, and frequencies uj\ i2 and uj 3 (thus left out is sum frequency 
generation and the generation of even higher harmonics). Then, the current density for held 
1 is given by 

Jx{wi,k u ) = a (1) (u;i, k lx )E x (ui,k lx ) + a (2) (u;i, k lx ] uj 2 , k 2x , u 3 , k 3x )E x (uj 2 , k 2x )E x (uj 3 , k 3x ). 

(S4) 

Here, E x (ui,ki x ) = t\E\j sin di is the total parallel held for i — 1 at the graphene layer, 
which we have written in terms of the incident held and transmission coefficient, cd 1 ) j s the 
linear conductivity function, while cd 2 )(u;i, ki x ;u; 2 , k 2x ,uj 3 , k 3x ) is the second-order nonlinear 
conductivity functions relating the current density generated at cui, ki x given helds at u 2l k 2x 



FIG. S7: Illustration of p-polarized electromagnetic fields propagating in the x-z plane and inter¬ 
acting with graphene. The fields i = 1,2 consist of incident, reflected, and transmitted components, 
with the directions of propagation and polarizations indicated by the red and black arrows, respec¬ 
tively. The angles of incidence and transmission are 8 and cj). The incident field is assumed to 
propagate in vacuum (refractive index n = 1), while the graphene sits on top of a substrate with 
index n (possibly frequency dependent). 






and ui 3 , k 3x . Similar expressions as Eq. (S4) can be written down for the current density 
at lO i,k ix (i = 2,3), and we use an analogous set of conventions to indicate the fields and 
conductivities at other frequencies and wavevectors. In what follows, we will also adopt 
the more compact notation = a^ 2 \u 1 , ki x ;u> 2 , k 2 x,u 3 , k 3x ), where the dependence on 

wavevectors and input frequencies is understood. 

The substitution of Eqs. (S3) and (S4) into Eq. (S2) (along with analogous equations 
for the other fields i = 2,3) yields a set of nonlinear equations relating the transmission 
coefficients and incident fields, 

ti = t[ L) 1 ~ t[ L) 4 L) aW (cai)o' (2) (u> 3 ) sin fa sin 2 <j) 2 sin <j) 3 , (S5) 

t 2 = 4 L) 1 _ ^ ElI \ t ( 2 L) i L) *^ 2 \co 2 )a^*(uj 3 ) sin 2 sin 0 2 sin* (fi 3 . (S 6 ) 

( 2 ceoj" 

Here t^ — r-. a ... — is the linear transmission coefficient, and the angle of 

n sin 6/+sin0+ (erf/ceo ) smt/sm0 7 ° 

the generated held is defined via k$ x = cos 03. The complex in-plane held amplitude 
generated at the difference frequency and at the position of the graphene layer z = 0 is given 

^ (L) 

E 3x = - EuE* 2I a (2) (u 3 ) sin 0i sin 0 2 sin 0 3 . (S7) 

zceo 

While Eqs. (S5) and (S 6 ) may appear somewhat complicated, here we note their main 
features. First, we note that the input held amplitudes, the beam angles, and the linear 

optical properties of the system are generally known. Thus, on one hand, given a theoret¬ 

ical model of the nonlinear conductivity these equations can be solved to obtain the 
predicted changes in transmission and reflection of the input beams, due to the generation 
of plasmons at CU 3 , k 3x . The plasmon held amplitude itself can be found from Eq. (S7). On 
the other hand, even absent a theoretical model, if changes in transmission or reflection of 
the incident helds are experimentally measured, one can attempt to invert these equations 
in order to obtain an experimentally inferred value of 

The description above generally holds regardless of the values of u 3 , k 3x . It is particularly 
interesting, however, to focus on the case where they align with the plasmon dispersion 
relation. In the small wavevector limit and for frequencies smaller than twice the Fermi 
frequency, u> < 2up, the linear conductivity is well-approximated by the Drude model [S5], 


cr(u) ~ 


ie 2 uip 
nhu + i'y 


(S8) 
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Here we have included a phenomenological damping term 7 . The plasmon dispersion relation 
can be found by solving for the pole of the linear transmission coefficient, t^ L \ For simplicity, 
we will momentarily consider the case of a substrate with frequency-independent refractive 
index, so that the role of plasmon damping can be more clearly identified. In the absence 
of losses, the dispersion relation is found to be 

, (i + 

™x,p , ? 

\OLCL)Jf 

where a ~ 1/137 is the fine-structure constant. In the presence of losses, choosing u> 3 , k 3x to 
lie on the plasmon dispersion relation yields a linear transmission amplitude of 
where Q = CJ 3/7 is the plasmon quality factor. Under these conditions, one thus sees from 
Eq. (S7) that the held intensity experiences a resonant enhancement of \E 3 \ 2 oc Q 2 . A similar 
resonant effect appears in the transmission and reflection coefficients of the incident fields. 

While the Eqs. (S5) and (S 6 ) can in principle be inverted to infer given experimental 
data for reflection or transmission coefficients, in the present experimental setup this proce¬ 
dure can only be done semi-quantitatively due to a number of unknowns. First, the signal 
lies significantly above the noise floor only near the plasmon dispersion relation, and only 
a limited number of beam angles are investigated. This makes it difficult to infer a specific 
wavevector and frequency dependence of the nonlinear conductivity (fundamentally, there 
must be a dependence on wavevector, as otherwise cd 2) = 0 for a centrosymmetric material). 
Furthermore, the experiment employs pulses whose bandwidths are significantly larger than 



the plasmon linewidth. Given that, here we aim to reach a conservative estimate for the 
strength of while we anticipate that future improved experiments (such as with longer 
pulses and nano-structures) and theoretical models will enable more detailed comparisons. 


The full conductivity function at zero temperature is given by [S5] 

. . ie 2 ojf e 2 |A. ^ i , u — 2ujp 

a M = ~ ■ + 7 T © ^ - 2 u F ) + - log ——-— 

7T/iU+«7 4/1 7T LO + ZUJf 


= -J —— + TT 0 (cn — 2u f ) + -log 
Tin u + ?7 An l 7r 


(S10) 


where Q(x) is the Heaviside step function. The Fermi energy hiO f ~ 0.5 eV of graphene 


is significantly lower than the pump and probe photon energies of ~ 3 eV. At these fre¬ 


quencies, the linear conductivity of graphene is nearly frequency independent and real, 
<7 ( 1 l(o;)/(ceo) ~ vra, which we use to obtain the pump and probe linear reflection coeffi¬ 
cients. Furthermore, we take the simplest possible function for the nonlinear conductivity, 
crl 2 l(o;) = i\o^ 2 \tjj 2 )\(u/ 112 ), where \<j ( ' 2 ' 1 (d 2 )\ is a single fitting parameter (the probe fre¬ 
quency CO 2 is fixed in the experiment). With this choice of function, graphene would be 
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equivalent to a nonlinear material with a frequency-independent bulk nonlinear suscepti¬ 
bility of x® = — ia^iu) / (ue^t) = |cr (o^ 2 ) I /(^ 2 ^ot), where t is the effective thickness of 
graphene. Inserting this nonlinear conductivity into Eqs. (S5), (S6), and (ST), we find that a 
value of |cr( 0 ^ 2 )| ~ 2.4 x 1CT 12 A-m/V 2 produces a good qualitative fit to the experimental 
data. 

We now discuss how to obtain the conversion efficiency of pump photons to plasmons. 
The number of photons dissipated per unit area and time by the field at the difference 
frequency consists of two terms, T d = T dg + T ds . The first term consists of damping 
from the graphene layer due to the real part of its conductivity, and is given by T dg = 
(Re cr^^oq)) \E 3x \ 2 /fihujs). The second term is due to damping from the substrate, due to 
the fact that at low frequencies its refractive index is complex. This contribution is given 
by T dtS = 4 fi | Im(fc ^ tan 0 3) | (Im n 2 (u 3 ))\E 3x \ 2 (l + | cot 1 2 )- On the other hand, the incident 
photon flux in the pump field is T in = / 1 sin 0 1 /(/ku 1 ), where I\ is the pump intensity. 



0 5 10 15 20 
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FIG. S8: The numerical solution for the pump conversion efficency, 77, as a function of wavevector 
and difference frequency. 


11 








Generally, the amplitude of the generated plasmon held rate will depend on both the pump 
and probe intensities. However, at the level of individual photons, the process is that an 
incoming pump photon gets converted to a plasmon (assisted by stimulated emission of a 
photon into the probe). Thus, we define the conversion efficiency relative to the pump alone. 
In steady state, the rates of photons dissipated and generated at the difference frequency are 
equal, and thus the overall conversion efficiency of pump photons to plasmons is r/ = 1^/1^. 
This efficiency is shown in fig. S8 for the range of frequencies and wavevectors used in our 
experiments. Using the estimated value of \a^(u 2 )\, we fold that the conversion efficiency 
for the experimental arrangement shown in fig. 3(b) at the point of maximum signal is 
approximately r) ~ 6 x 1CT 6 . 
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